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Abstract 

In this paper we compare three Perfectly Matched Layer (PML) treatments by means 
of a series of numerical experiments, using common numerical algorithms, computational 
grids, and code implementations. These comparisons are with the Linearized Euler Equa- 
tions, for base uniform base flow. We see that there are two very good PML candidates, 
and that can both control the introduced error. Furthermore, we also show that corners 
can be handled with essentially no increase in the introduced error, and that with a good 
PML, the outer boundary is the most significant source of error. 


1 Introduction 

Aeroacoustics research is conducted with the combined use of theoretical analysis, experimental observations, 
and computation. Computational AeroAcoustics (CAA) is growing as a result of increases in the capability 
of numerical methods and computer hardware, and is increasingly valuable both as a source of insight and 
as a tool for design [3, 17, 22, 27]. Artificial boundary treatments are a critical element for CAA because 
they permit the restriction of a numerical domain while causing an error in the numerical solution that 
is intended to be limited. Various reviews of the extensive work done to date on this problem have been 
published, including [2, 5, 6, 12, 14, 15, 23, 26]. 

Hagstrom [13, 18, 4] has analytically developed a PML treatment for artificial radiation boundaries in 
the case of the LEE with uniform base flow. A widespread and general alternative approach is to use a 
sponge zone, as in [1, 19, 20, 21]. Kami [16] discusses possible forms for damping in a reflectionless sponge 
zone for a first order hyperbolic system. In this paper we use two PML sponge layer treatments. The 
two sponge layer treatments can be extended to the LEE with nonuniform steady and unsteady base flows, 
and to the Compressible Navier Stokes Equations. These three PML treatments are compared by means 
of a series of numerical experiments, using common numerical algorithms, computational grids, and code 
implementations. Two further series of experiments with one of the sponge layer treatments show how to 
treat corners where two PMLs overlap, and that the outer boundary treatment can cause large errors even 
with a good PML treatment. The numerical experiments are designed to isolate three aspects that are 
important for the actual use of the PML: (1) a simple PML at a single boundary, with a layer that is wide 
enough so that the propagated disturbance does not reach the outer boundary; (2) a combination of two 
PMLs at a common corner, with a layer that is wide enough so that the propagated disturbance does not 
reach the outer boundary; and (3) a PML combined with a treatment for the outer boundary of the layer. 
Within the context of these examples, we see that there are two very good PML candidates that can both 
control the introduced error, that corners can be handled with essentially no increase in the introduced error, 
and that with a good PML, the outer boundary is the most significant source of error. 
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2 PML for the LEE in two space dimensions 


The Linearized Euler Equations (LEE) are obtained by linearizing the Euler Equations about a base solu- 
tion. We will consider two dimensional cases where the base solution {pb,Ub,Vb, Pb} is steady. We use a 
nondimensional form of the governing equations, and note that a steady parallel flow with constant pressure 
and density is a solution for the Euler Equations. 

The simplest case for the LEE in 2D is with a constant base solution. We express the constant base 
velocity components are expressed in terms of Mach number, ( , Vb) = ( M x , M y ). The perturbation density 
is uncoupled from the perturbation pressure and velocity, and will be neglected. In this case, the nondimen- 
sional form of the LEE for perturbation pressure and velocity components can be written as 


dp , , dp , , dp du dv 
- + M X — + M,- + — +- = 0, 

(i) 

du , , du , , du dp 

at + M ‘di + M, di + di = °’ 

(2) 

5 + mA + m.A + I^o. 

dt dx v dy dy 

(3) 

If we let the unsteady perturbation quantities be 
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(4) 


then this simplest form of the LEE in two space dimensions can be written as 


dU 

~dt 


where A and B are the constant matrices 
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( 6 ) 


A Perfectly Matched Layer (PML) boundary treatment is designed for this system by ensuring that the 
plane wave solutions for the LEE in the numerical domain perfectly match the plane wave solutions for the 
adjacent PML domain along the interface between them. 

We shall consider the LEE with a pressure source term, with the modified pressure equation 


dp 


dp 


dp du dv 


^7 + M,— + M y — + — + — = asin[27 rut] exp [-(fi x x 2 + (3 v y )} = g[x , y, f), 


dt 


dy dx dy 


(7) 


for all ( x , y) in the numerical domain. With this source term, the LEE can be written as 


du du du 

~dt +A ~d^ +B ~d^- G ' 


(8) 


where G = ( g , 0, 0) T inside the numerical domain. For the numerical experiments reported in this paper, we 
always take 

a =0.01, (3 X = 36 = f3 y , v=l. (9) 

Note that exp[— 36 * 2 2 ] ss 10~ 63 . This radiating source produces a propagating disturbance that is O[10 -4 ] 
in the numerical domains that we use. We do not include the source term in the pressure equation inside the 
PML, and take a large enough numerical domain so that the pressure source and all of its relevant derivatives 
are effectively 0 at the interface. 
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2.1 The Hagstrom PML 

Hagstrom [13] has formulated a PML boundary treatment for this system, which has been further investigated 
by Motamed [18] and Froling [4]. For an artificial boundary layer that is perpendicular to the x axis, this 
PML treatment can be expressed as the system 


dtJ 

~dt 



pa(x)U 


dU 

B— + W = 0, 
dy 


(10) 


dW 


M 


dW dU - \ 

— 1- cr(x)W + a(x)A [ — — b pa(x)U ] = 0, 

dy \ ox I 


( 11 ) 


where W = (w p ,w u ,w v ) T is a vector of auxiliary variables, a(x ) is a damping profile, and 


9 = 


M x 

i -Mr 


(12) 


This PML system is used in a layer that is joined or matched to the computational domain of interest. If the 
numerical domain is (x,y) € [Xi,X r ] x [Y],, Y t ], then a PML domain on the right hand boundary would be 
for {x, y) £ [X r , X r + D r \ x [Yj, Yj], where D r is the layer width. In such a layer, a typical damping profile 
from [4] is 

" w = “( £ ^r 1 ) , (is) 

where a is a constant, and usually n = 3, 4, or 5. We shall use 


( x ) = c A 

( ) c {i + {x-x r yy' 


and will call the Hagstrom PML with this damping profile PML A. 


(14) 


2.2 The Kami directional damping PML 

Kami [16] has given a very nice eigensystem analysis of PML treatments with filtering operators for both 
slowing down and damping outgoing waves. By considering diagonalization in a general direction, Kami 
gives one particular form for directional damping in an artificial boundary layer as the system 

+ A~?r- + + 5(cos(9)A + sin(9)B)U = 0, (15) 

at ox ay 

where 9 is the angle of the damping direction, and <5 is the damping rate. For a damping PML perpendicular 
to the x axis, this reduces to the system 


dU BU dU 

~j 7 ~ + A — — b B— zfc SAU — 0, 
dt ox dy 


(16) 


where +<5 is taken on the right, and — <5 is taken on the left. This PML system with <5 = a will be called 
PML B. We note that the damping terms are 


d p = ±cr(M x p + u), d u = ±cr(M x u + p), and d v = ±crM x v, (17) 

in the equations for p, u, and v, respectively, where + is taken on the right, and — on the left. Notice that 
if M x = 0, then the damping terms for p and u are proportional to each other, and that v is undamped. 
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2.3 A simple PML by analogy 

The simplest example of a decaying dynamic is radioactive decay, governed by the ordinary differential 
equation 

f- + cf = °, (18) 

with c > 0. If we consider the LEE as an evolutionary operator that is analogous to a time derivative, then 
by analogy a simple damping system is 


dU dU dU - n 

-q 7 + A— h B— h oU — 0. 

at ox ay 


(19) 


We shall call this PML treatment PML C. We note that the damping terms are 


d p = ap, d u = au, and d v = av, 


(20) 


in the equations for p, 
underlying steady flow, 
case M x = 0. If 


u, and v, respectively. Notice that these damping terms are independent from the 
and that the p and u terms are the opposite of the damping terms in PML B for the 
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(21) 


then PML C is obtained from the characteristic analysis of a general damping form. What we call PML C 
has been discussed in [20], where the damping term aU is called a relaxation source term, and the damping 
coefficient a — 1/r, where r is a relaxation time. 


2.4 Characteristic based PML 

Kami [16] observes that the damping rate need not be the same for the different eigenfunctions of the 
hyperbolic system, but that the damping terms should not alter the eigenvectors of the hyperbolic system 
with respect to the damping direction. A characteristic analysis similar to Thompson’s [24, 25] can be used to 
obtain a general damping layer. The characteristic variables or Riemann invariants for the LEE diagonalized 
in the x direction are 

R = p + u, L = p — u, and v, (22) 

with propagation speeds M x + 1, M x — 1, and M x , respectively. If 0 < M x < 1, then R is right going, L is 
left going, and v is not left going. The LEE with damping terms can be expressed in terms of these variables 
as 

OR OR OR dv 

— + ( M x + l)— + M y — + - + a R ( M x + l)R-0, 

dL _ dL , , dL dv 

at + (M ' - ' 


+ (M x - 1 )— + M y gy ■ dy 


+ v; — b o’ l{M x — 1 )L — 0, 


dv dv dv 1 dR 1 dL 
_ + M x — + M y — + + <r v M x v = 0, 

where cr^j, a l, and a v need not be the same. The damping terms for R, L, and v are 

dn = o r (M x + 1 )R, d L = o l (M x - 1 )L, and d v = o v M x v, 


(23) 

(24) 

(25) 

(26) 


respectively. If this form of the equations is recast in terms of the primitive variables, then the damping 
terms in the primitive variable equations are 
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and 

d v = a v M x v = cv, 

respectively. The LEE with these primative variable damping terms can be written as 

dU A dU BU n 

^ 1 " DU — 0 , 

at ox ay 


with the symmetric damping matrix D is 


aR - aL M x 


£R_ l£L_' 
2 


D= ( + (£a±£tM x + 2^) 

0 0 


0 ^ 


( a 

6 

0 \ 

0 


b 

a 

° 

a v M x / 


K o 

0 

c 


If we consider the primitive damping coefficients a, 6, and c to be independant, then 
a + b = a R (M x + 1), a—b= a l(M x — 1), and c = a v M x . 


(29) 


(30) 


(31) 


(32) 


The eigenvectors for A are: E R = (1, 1, 0) T with eigenvalue A ar = M x + 1; El = (1, —1, 0) T with eigenvalue 

A al = M x — 1; and E v = (0, 0, 1) T with eigenvalue A a v = Mx- The eigenvectors for D are the same as for 

A , with eigenvalues 

A dr = a + b = a R (M x + 1) = c t r \ar , (33) 

A dl = a - b = <jl{M x - 1) = ct^Aal, (34) 

and 


A Dv — C — &vM x — <J v \avi 


(35) 


respectively. 

Both PML B and PML C use just one damping profile, with cr^ = <jl = <x v = If we set a = M x a = c, 
and b = a, then PML B is recovered, with 


f Mx l 0 A 

D b = a(x) 1 M x 0 (36) 

V 0 0 M x ) 


If tr > 0 in a PML B, then R is damped as it goes right, L is amplified as it goes left, and v is either damped 
if M x > 0, or unchanged if M x = 0. Note that if M x = 0, then the damping matrix for PML B becomes 


D b \m x =o = o-(x) 


0 1 0 \ 
10 0 
0 0 0 / 


If a = a = c and 6 = 0, then PML C is recovered, with 


D c = a( x) 


1 0 0 \ 
0 10 
0 0 1 / 


(37) 


(38) 


If tr > 0 in a PML C, then i?, L and v are all damped at the same rate. For a PML perpendicular to the y 
coordinate axis, we note that 

/ d 0 e A 

£= 0 / 0 (39) 

\ e 0 d ) 

is an acceptible general form for a damping matrix similar to D for the case in the x direction. 
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3 Numerical Algorithm 

The Hermite/Cauchy-Kovalevsky/Taylor method (see [7], [8], [9], [10], [11]) is used for all of the experiments 
that are reported in this paper. This method uses two staggered uniform grids, offset by half a grid spacing 
in both x and y. The first or base grid is used for initial data, and the numerical solution at each full time 
step. The second or offset grid is used for the numerical solution at each half time step. The time evolution 
through a half time step from one grid to the next starts with multidimensional Hermite spatial interpolation 
of the local data on the current grid stencil, with a tensor interpolant for a rectangular grid, it then uses 
the governing equations for Cauchy-Kovalevsky recursion to produce all of the required time derivatives in 
space and time from the spatial derivatives produced by the interpolant, and it ends by propagating through 
a half time step with a Taylor series. This approach is not simply an algorithm, but is actually a method 
for developing or specifying numerical algorithms. This approach has had various names, but in this paper, 
we will call it the Hermite/Cauchy-Kovalevsky/Taylor method, or the HCKT method. 

The spatial interpolant with data on one grid is expanded about the center of the stencil, which is a 
point in the alternate grid, so that the spatial derivatives from the interpolant with known data on one grid 
are used by the Cauchy-Kovalevsky recursion to compute time derivatives at a grid point of the alternate 
grid. A simple method is used for the numerical experiments that are reported in this paper. This method 
uses a 2 x 2 local stencil, and at each of the four grid points in the stencil, for each of the dependant 
variables, we carry data for the variable itself, and its x, y , and xy derivatives. The spatial interpolant is 
a third order Hermite tensor interpolant on a stencil that is two points wide, and it simultaneously and 
consistently estimates all the derivatives of order {( in , n) : 0 < m,n < 3}, with up to the sixth order (3, 3) 
cross derivative. Other methods could be specified with different stencils and different data, but we use only 
one simple method since the PML issues are the focus of interest here. 

The Cauchy-Kovalevsky recursion is used with the local interpolated data surface considered to be initial 
data, and it can produce time derivative data of any order. The recursion formula comes directly from the 
governing equations, and in the simple constant coefficient case, it is 

gi+j+k+1 -fj Qi+j+k+ljj Qi+j+k+lfj 

dx'dyidt k+1 dx l+1 dyidt k dx' l dy^ +1 dt k ' ^ ^ 

for k > 0, and for all required i and j, depending on the interpolant and the upper range of k. The time 
derivatives that are obtained from this recursion are used in a simple Taylor series time advancement at the 
stencil center, so that known data on one grid produces time derivaties for a solution on the alternate grid. 
For the LEE with constant coefficients, and with the method that we use that has a third order interpolant, 
if the recursion is carried through the sixth order in time, and if all of the time derivatives are used for 
the Taylor series time advancement, then the locally interpolated two dimensional data surface would be 
propagated in time exactly. Since the propagation of the local two dimensional data surface is exact, this 
method correctly captures propagation along characteristic surfaces in multiple spatial dimensions, and in 
this sense, it can be said to be a correct generalization of the method of characteristics to higher dimensions. 
For the results presented here, we only compute and use time derivatives through the fourth order for the 
third order interpolant, so that the propagation is not exact, but it still is correctly along the characteristic 
surfaces to the fourth order in time. The numerical method is third order in time and space for p, u, and v, 
it is second order for their x and y derivatives, and it is first order for their xy cross derivative. Overall, this 
method is third order in time and space for the physical variables, with diffusive leading order error terms. 

By changing the stencil width, or the depth of data at each stencil point, various other interpolants can 
be used, with arbitrary order. By changing the recursion to a different set of governing equations, a different 
physical system can be simulated. Any system of partial differential equations can be approximated with 
this HCKT algorithm development method as long as the governing equations can be written in the form 

dV 

^ = p (V), (41) 

where V is the vector of dependant variables, and where F is any partial differential operator defined in 
terms of V and its spatial partial derivatives. For example, the various PML equation systems specify 
governing equations in the layers that are different from the LEE that govern in the numerical domain. 
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We have used the HCKT method to develop algorithms for a variety of governing equations, including the 
advection equation, the heat equation, the convection diffusion equation, a first order dispersive system, 
the LEE, Burger’s equation, and the Compressible Navier-Stokes Equatiob, and these algorithms have had 
accuracy up to the twenty first order in both space and time. Algorithms can be specified by the HCKT 
method by simply replacing the recursion equations or subroutine. A comparison of accuracy and efficiency 
for algorithms up to the fifteenth order for a simple one dimensional problem is given in [10]. 


4 Experiments With One PML Boundary 

In this section we presents results from a series of numerical experiments using PML A, B, and C. The 
computations that we present here are all on a rectangular numerical domain, 


n N = {{x,y)e [X h X r \x [Y b ,Y t ]}, 
with an additional PML layer on the right 

n R = {(*, y) e [X r , X r + D r ] x [Y b , Y t ]}, 


(42) 


(43) 


where D r is the layer width. All of these computations are for the LEE with a Gaussian pressure source. 
The governing equations in fijv are 


where 


dU 

~dt 


A 


f + B^. a . fa .0,0r, 


g = 0.01sin(27rf)exp(— 36(a; 2 +y 2 )). 


(44) 

(45) 


The source is negligible except near the origin, it radiates a signal away from the origin, and it is not included 
in the PML equations. Since we have an interest in acoustic applications, we consider the propagation 
medium to be scaled by the standard atmosphere, with the pressure of the base flow to be one atmosphere, 
and the amplitude of the pressure disturbance has been chosen to produce a radiating signal that is O[10 -4 ]. 
The frequency is v = 1, so that the wavelength of the radiated sound is one spatial unit. If we consider 
the unit time to represent Is and the frequency to represent 1Hz, then the spatial unit corresponds to 
approximately 340m, while if we consider the spatial unit to correspond to lm, then the time unit scales 
so that v = 1 corresponds to a frequency of approximately 340Hz. The numerical domain is always chosen 
large enough, and the elapsed simulation time short enough, so that in these particular experiments, the 
radiated signal only leaves the numerical domain on the right at the interface between the numerical domain 
Hat and the PML layer Qr. In other words, we use Big Enough Domain (BED) boundary conditions along 
the numerical domain limits at X), Y b and Y t . The damping profile that we use is 


a(x) 


(x - X R ) 6 

C (l+(cr-cr«) 2 )3’ 


(46) 


for X r < x < X r + D r , where D r is chosen large enough so that we have a BED boundary condition at 
X r + D r , and no signal exits the PML layer CIr at its outer boundary. Note that we are now using c as a 
real constant, and not as a function in a damping matrix, as earlier. A third order algorithm is used, and 
the damping profile has been chosen first of all so that a and sufficient of its derivatives are 0 at the interface 
between the numerical domain and the PML layer to ensure a completely smooth seam along the interface, 
through the order of the leading order truncation error terms of the numerical algorithm. This means that 
one sided characteristic matching at the interface is unneccesary, since there is no difference up to the order 
of the numerical method between the LEE and the PML system at the interface. We simply change the 
recursion routine from to f }r, and we use the numerical domain recursion on the interface at the half 
time steps when a cell center is at Xr. The damping profile has also been chosen to smoothly vary from 0 
along the interface to an essentially constant value deeper inside the PML layer. All computations are with 

A*=1 = A y, At=l (47) 
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In all cases, the reported data has been obtained by comparison with a reference solution computed with 
c = 0, so that the PML layer Hr is actually just an extension of the numerical domain Qjv, with a BED 
boundary condition at the outer boundary of the layer. We have always taken X r = 2, with data obtained 
as 

E L {t) =ma,x{\p L (x,yJ) - p r ef{x, yj)\ : {x,y) GH n } (48) 

where p re f is the pressure from the reference solution, and p l is the solution using PML L , for L = A, 
B , or C . The reported data is not an error in the sense of a comparison of a numerical estimate and an 
exact solution, but it is a measure of the difference betweeen two numerical estimates, one on a Big Enough 
Domain, and one with a PML. The reported data indicates only the effect on a numerical solution that is 
produced by truncating the numerical domain with a particular PML. 


Table 1: Maximum absolute disturbance from one PML on the right. 
PML A, B, and C with Az = 1/24 = Ay, At = 1/48. 

Data with M x = 0.0, 0.4, or 0.8 at t = 10. 


c 

Ea{ 10) 

E b ( 10) 

E c { 10) 


M x = 0.0 

0.001 

5.7686 x 10 _i4 

7.3605 x 10 _i2 

6.4692 x 10" y 

0.010 

5.7710 x 10" 13 

7.3643 x 10 _ii 

6.4552 x 10" 8 

0.100 

5.7946 x 10" 12 

7.4020 x 10" 1U 

6.3175 x 10" Y 

1.000 

6.0052 x 10" 11 

7.7960 x 10" y 

5.2127 x 10" y 

10.00 

6.7827 x 10" 1U 

1.5901 x 10"' 

2.7457 x 10 -t> 


£ 

II 

o 

0.001 

7.1998 x 10" 13 

1.4915 x 10" 12 

2.9545 x 10" y 

0.010 

7.2031 x 10" 12 

1.4934 x 10" 11 

2.9506 x 10" 8 

0.100 

7.2357 x 10 _il 

1.5122 x 10 _1U 

2.9115 x 10" 7 

1.000 

7.5906 x 10 _iU 

1.7337 x 10" y 

2.5720 x 10" y 

10.00 

1.0285 x 10" 8 

7.9125 x 10" 8 

1.4206 x 10 -t> 


M x = 0.8 

0.001 

1.3805 x 10" 13 

6.0922 x lO" 14 

3.3349 x 10 _iU 

0.010 

1.3848 x 10" 12 

6.1014 x 10 -13 

3.3219 x 10" y 

0.100 

1.4283 x 10" 11 

6.1966 x 10" 12 

3.2596 x 10" 8 

1.000 

1.8633 x 10" 1U 

7.1884 x 10" 11 

2.9821 x 10" 7 

10.00 

Blew Up 

3.0027 x 10" y 

1.9087 x 10 _y 


The computations reported in Table 1 are for three cases of uniform mean flow, with M x = 0.0, 0.4, or 
0.8. In the first case with M x = 0, the evolution of p is also governed by the standard wave equation. In all 
cases the reported data is at t = 10. The numerical domains and PML layers are 

H N = [-11, 2] x [-11, 11], and H r = [2, 11] x [-11, 11], for M x = 0.0, (49) 

H n = [-7,3] x [-11,11], and Hr = [3, 15] x [—11, 11], for M x = 0.4, (50) 

H n = [-3,4] x [-11,11], and Hr = [4,19] x [-11,11], for M x = 0.8. (51) 

In order to scale the data in Table 1, the pressure solutions with the expanded BED numerical domains are 
O[10 -4 ] in all cases. Notice that the disturbance of the numerical solution that is produced by introducing 
a PML ranges in absolute terms from O[10~ 14 ] to 0[1O -6 ], or in relative terms from O[10~ 10 ] to 0[1O -2 ]. 
None of the disturbances in the numerical domain solution can be seen at the scale of the solution if they 
are relatively smaller than O[10~ 3 ]. We note the evident superiority of PML A and B when compared to 
the simple PML C. We also note that PML A works better at lower convection speeds, while PML B works 
better at higher convection speeds. From this series of comparisons it appears on the whole that PML A 
and PML B are reasonably comparable, with the advantage of PML A at lower speeds being somewhat 
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better than the advantage of PML B at higher speeds. Recall that PML A uses a set of auxiliary variables 
and equations in the PML, while PML B uses only the primitive variables, so that PML A requires twice 
the amount of storage and has twice as many equations in the PML layer. Also recall that PML A is 
specifically taylored for the LEE with constant coefficients, or for uniform flow, while PML B can be readily 
extended to nonuniform or even unsteady flows by using variable coefficients in the PML B system. The 
numerical domain solution disturbances scale with the damping coefficient c, and the effect of the PML can 
be made essentially as small as desired with a small enough damping coefficient. This shows that these 
PML treatments actually do smoothly match across the interface, and that their effect on the solution in 
the numerical domain is controllable. We must recall that using a smaller damping coefficient to control the 
numerical domain disturbance will require a wider layer, so that getting a more accurate solution has a cost, 
as is only to be expected. 


5 Treating a Corner Between Two PML Boundaries 

Treating a corner between two PML boundaries has been a distinct problem with its own difficulties. We 
note again that Kami [16] has given a particular form for directional damping in an artificial boundary layer 
as the system 


+ <5(cos(0)A + sin (0)B)U = 0, 

at ox ay 


(52) 


where 0 is the direction of damping, and <5 is the damping rate. Notice here that <5 is the same damping 
rate for both AU and BU. For damping in the direction of a coordinate axis, this directional form reduces 
to PML B in that direction. We consider the corner between a PML B in x on the right, with domain 
{ A R = [X R , X R + D r ] x [Y b ,Y t ] and equations 


dU nu BU n 

~dt +A ^ + B ^ + (TR{x)AU - Q ' 

and a PML B in y on the top, with domain fly = [Xr, Xr\ x [Yt, Yr + Dr] and equations 

f+4+4+-<^ 


(53) 


(54) 


Note that <Jr{x) and <JR{y) need not have the same functional form. In this upper right corner we use a 
blended PML B form, with domain VLrr = [X Rl Xr + Dr] x [Yr, Yr + Dr] and equations 


dU_ A dU_ + R dU_ (x - x r)(7 r(x) 

dt dx dy ^{x - x R ) 2 + (y - y T ) 2 


AU - 


(y-yT)cr T (y) 


\/(x - Xr) 2 + (y - y T ) 2 


BU = 0. 


(55) 


This blended PML B form in U,rr is not a directional damping except along the bottom and left edges of 
the corner PML domain, even though 


cos (9(x, y)) = 


(x - Xr) 


\/{x- Xr) 2 + (y- y T )‘ 


= , sin (0(x,y)) = 


(y - Vt) 


^(x-x R ) 2 + {y- y T ) 2 ' 


(56) 


where 9(x , y) is the angle to the point ( x , y) relative to the U,rr corner point ( Xr , Yr). Many PML treatments 
can be succsefully blended like this at a corner, for example, any of the PML with the general characteristic 
form dicussed above, and in particular PML C. It is not immediately clear how to implement a blended 
corner treatment like this for PML A with its auxiliary functions. 

The computations reported in Table 2 are all for a uniform mean flow, with columns for each of the 
cases ( M x ,M y ) = (0.0, 0.0), (0.4, 0.0), (0.8, 0.0), and (0.8, 0.4). In the last case, the convection speed is 
| (M x , M y ) | « 0.89, and is still subsonic. The numerical domain is 


{l N = [Xi,X r \ x [Y b ,Y t ], 


( 57 ) 
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Table 2: Maximum absolute disturbance at t = 10 from a top right blended corner. 
PML B on the right and the top with Ax = 1/24 = Ay, At = 1/48. 


c 

(M x ,M y ) = (0,0) 

(M x ,M y ) = (.4,0) 

(M x , M y ) = (.8,0) 

(M x ,M y ) = (.8, .4) 

0.001 

7.3605 x 10 -12 

6.5968 x 10 -i:,i 

4.7396 -12 

9.6336” 13 

0.010 

7.3643 x 10 _il 

6.5988 x lO" 11 

4.7413 -11 

9.6505-^ 1 

0.100 

7.4020 x 10" 1U 

6.6184 x 10" 10 

4.7580" 1U 

9.8451' 11 

1.000 

7.7960 x 10” y 

6.8662 x 10" y 

4.9670 _y 

1.2351 _y 

10.00 

1.5901 x 10" 7 

1.5922 x 10" 7 

1.4395 -7 

7.5069- 8 


the PML layers on the right and top are 

n R = [X r , X r + D r ] x [Y b , Yt], fl T = [Xi, X r ) x [Y t , Y t + D t ], (58) 

and the blended corner PML domain is 

n RT = [Xr, X r + Dr] X [Y b , Y t + D t \, (59) 

where, for ( M x ,M y ) = (0.0, 0.0) 

X t = -11, =2, D r = 9, Y b = -11, Y t = 2, and D t = 9, (60) 

for (M x , M y ) = (0.4, 0.0) 

Xi = -7, X r = 3, D r = 12, Y b = -11, Y t = 2, and D t = 9, (61) 

for (M x , M y ) = (0.8, 0.0) 

Xi = -3, = 4, D r = 15, Y b = -11, Y t = 2, and D t = 9, (62) 

and for (M x ,M y ) = (0.8, 0.4) 

Xi = -3, = 4, D r = 15, Y b = -7, Y t = 3, and D t = 12. (63) 


In all cases both <j r (x) and (Jt(v) have the same functional form that was used above for a single PML 
layer on the right. As for Table 1, the reported data is the maximum absolute difference in the numerical 
domain at t = 10 between the solution with the PML treatments and a solution on a large domain with 
BED boundaries, so that the propagated disturbance in the BED numerical domain never reaches the outer 
boundaries. The pressure solutions with the expanded BED numerical domains for all of these blended corner 
cases are O[10 -4 ]. In all cases, the results in Table 2 with PML B treatments on the right and the top, and 
a blended PML B corner treatment are very similar to the results in Table 1 where there is only one PML B 
on the right. Notice that the data for the ( M x , M y ) = (0.0, 0.0) case in Table 2 is identical to the comparable 
data in Table 1 for M x = 0.0. The disturbances introduced with a blended corner PML treatment for the 
(M x ,M y ) = (0.4, 0.0) and (. M x ,M y ) = (0.8, 0.0) cases in Table 2 are of the same order as the disturbances 
for the comparable M x = 0.0 cases in Table 1, which would be the expected disturbance levels introduced 
by the single PML B treatment along the top of the numerical domain with M y = 0.0. The disturbances 
introduced with the blended corner PML treatment for the skew (M x ,M y ) = (0.8, 0.4) case in Table 2 are 
smaller than the comparable case in Table 1 for either the M x = 0.4 or the M x = 0.8 cases. We note that 
the disturbance introduced into the numerical domain solution by a single PML B layer is actually smaller 
as the convection speed toward the layer increases, and that the speed in this skew case is approximately 
0.89. The data in Table 2 shows that the blended corner PML B treatment introduces essentially no new 
disturbance into the numerical domain solution beyond what is caused by the one dimensional PML B layers 
on the right or the top of the numerical domain. 
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6 The Effect of the Outer Boundary 

Reflection from an outer boundary back through a PML into the numerical domain is another distinct 
problem with its own difficulties. In essence, the outer boundary of a PML layer presents the same problem 
as the original artificial boundary, only removed a distance from the numerical domain. What is a good 
transparent boundary condition anyway? The entire PML enterprise has been created because there is no 
easy answer to this question. The hope in using a PML is that a propagating disturbance will vanish before 
the outer boundary is reached. The problem in using a PML is that even though there is damping through 
the layer from the numerical domain towards the outer boundary, there can also be amplification from the 
outer boundary back towards the numerical domain. 

We consider here a range of three simple Boundary Conditions for an outer boundary for a PML. The 
first is trivial, simply setting all disturbance data to 0 at the outer boundary. We will call this first outer 
boundary treatment OBi . The second outer boundary treatment is to extrapolate R = p + u and v at the 
outer outflow boundary, and to set L = p — u = 0. If the outer boundary is at X ob = X r + D r , then the 
data at X ob _ i = X ob — Ax and X ob - 2 = X ob — 2Ax is interpolated with a cubic Hermite interpolant. This 
interpolant is used to extrapolate the data at X ob , with 


f{X ob ) = — 4/(X 0 (,_ 1 ) + 5 f(X ob _ 2 ) + Ax [4/ x (X o6 _!) + 2f x (X ob -i)\ , (64) 

and 

f x (X ob ) = 12 [~f(X ob _ r) + f(X ob _ 2 )] /(Ax) + 8f x (X ob - 1 ) + 5 f x (X ob _ 2 ). (65) 

Two interpolation problems are solved for each of f = R and / = v, one that gives / and f x , and another 
that gives f y and f xy . The physical variables p and u , and their x, y and xy derivatives, are all recovered 
at the outer boundary from the extrapolated data for R and the zeroed data for L. We will call this second 
outer boundary treatment OB 2 . The third outer boundary treatment is to propagate R and v with interior 
differencing, with L = 0. The propagation equations with damping are 


OR 

~dt 


,dR 

dx 


OR 


— + (M x + 1) — + My — + — + a R (M x + 1 )R - 0, 


dv 

dy 


( 66 ) 


L = 0, 


(67) 


dv , , dv , T dv IdR 

— + M x — — b My— — |- — — — h a v M x v — 0. 
dt dx dy 2 dy 


(68) 


Note that at this outer boundary the term (l/2)dL / dy has been dropped from the equation for v. We will 
call this third outer boundary treatment OB 3 . Note that the constraint L = 0 is stronger than the constraint 
L x = 0. 


6.1 PML B on One Side 

In this section we consider simulations with PML B on one side in order to see what can happen to a 
numerical solution because of the reflection of a signal from the outer boundary of a PML. We consider 
only the case ( M x ,M y ) = (0,0) with no convection. Each computation uses PML B, with outer boundary 
treatments OBi , for i = 1, 2, or 3. All computations have the PML interface at X r = 2. All of the domain 
configurations are designed so that the propagating signal does not reach the outer boundaries on the left, 
top, and bottom. 

The data in Table 3 is from three sets of computations, at two different times with two different PML 
widths. As before, the reported data is the maximum disturbance in the pressure in the numerical domain as 
calculated by taking the difference between the numerical solution with the PML and a numerical solution 
with a Big Enough Domain that has no layer, and consequently no interaction with any outer boundary. 
We note that the first and second sets of data have the same layer width at different times, and the second 
and third set have different layer widths at the same time. The first set is with PML width D r = 5 at 
t = 15, with Xi = 16 = Y b = Y t . The outer boundary of the PML is at X ob = 7, so that by t = 15 the 
signal will propagate to the outer boundary, and any reflection will have a chance to return back through 
the PML and reach as far as x = —1, well within the numerical domain. We note that all of the accuracy 
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Table 3: Maximum absolute disturbance from an outer boundary. 
PML B with Az = 1/24 = Ay, At = 1/48. 


c 

OB 1 

ob 2 


ob 3 



t = 15 and w = 

5 


0.100 

9.3051 x 10 _B 

8.3402 x 10 _B 


9.3046 x 10 _B 

1.000 

9.2273 x 10" y 

8.2684 x 10 -y 


9.2267 x 10 -y 

10.00 

4.7438 x 10 _B 

4.2690 x 10“ B 


4.7431 x 10 _B 



t = 25 and w = 

5 


0.100 

2.5062 x 10“ 5 

2.3514 x 10 -t> 


2.5062 x 10 -t> 

1.000 

2.4807 x 10" b 

2.3275 x 10" b 


2.4807 x 10" b 

10.00 

1.0290 x 10" b 

9.6544 x 10 -y 


1.0289 x 10 -t> 



t = 25 and w = 

10 


0.100 

4.0624 x 10 -y 

3.6134 x 10” y 


4.0623 x 10" 6 

1.000 

3.9743 x 10" y 

3.5361 x 10 _B 


3.9741 x 10" y 

10.00 

1.7670 x 10- b 

1.7670 x 10- b 


1.7670 x 10- b 


that was exhibited with this same PML B in Table 1, has disappeared, that all of the disturbances are now 
O[10~ 6 ]. It is worth noting that the data with c = 10 shows a disturbance that is half what is shown for the 
smaller damping coefficients. For the cases in Table 1, where the propagating signal never reached an outer 
boundary, the disturbances introduced by the PML scaled with the damping coefficients, so that the worst 
performance was with c = 10, and the effect of introducing a PML could be controlled and made as small 
as desired by choosing a sufficiently small value for c. The second set of data is with PML width D r = 5 
at t = 25, with Xi =27 = Y& = Y t . The outer boundary of the PML in this case is also at X 0 b = 7, so 
that by t = 25 the signal will propagate to the outer boundary, and any reflection will have a chance to 
return back through the PML and reach as far as x = —11. The distortions introduced by the PML in this 
set are approximately twice as large as in the first set above with the same layer width at t = 15. This is 
approximately linear growth in time. We note that the OB 2 extrapolation BC is slightly better in this case 
than the other two. The third set is with PML width D r = 10 at t = 25, with X/ = 27 = Yy, = Y t . The outer 
boundary of the PML in this case is at X 0 b = 12, so that by t = 25 the signal will propagate to the outer 
boundary, and any reflection will have a chance to return back through the PML and reach as far as x = — 1. 
The distortions introduced by the PML in this case are approximately an order of magnitude smaller than 
in the second case above with a layer that is half as wide at t = 25. This is an approximately linear decrease 
with layer width, and shows the advantage of increasing PML width. 

Considering all of the data in Table 3, it appears that all of the simple outer boundary treatments that 
we have tried are equally bad. The disturbances introduced by these simple outer boundary treatments 
completely obliterate all of the quite good accuracy that was obtained by the PML itself, as demonstrated 
by the data in Table 1, where an outer boundary was not significant. The disturbances are all in the absolute 
range of 0[1CC 6 ] to 0[1O -5 ], which is approximatrely 1% to 10% of the numerical solution with no boundary. 
In this regard, it is worth noting that an additional series of computations with PML width D r = 5 at t = 10 
produced disturbances for each value of c that were identical to the data in Table 1, irrespective of the 
outer boundary treatment, since the reflected signal did not have enough time to propagate back into the 
numericl domain. We see clearly that the outer boundary can be a significant source of contamination for 
numerical solutions with artificial boundary treatments. This is at least true for the three outer boundary 
treatments that we have tried here, which included the extremely crude OBi that simply sets all data to 0 
at the outer bounadry, the extrapolation BC OB 2l and the OB 3 that uses interior differenced propagation 
for the outgoing characteristic variables with the strong Thompson like condition L = 0 for the incoming 
characteristic variable. The results from these three outer boundary conditions were very close to each other, 
typically within about 10% of each other. The underlying issue undoubtedly is that the PML B both damps 
going to the right, and amplifies going to the left. However, it remains to be seen if a combination of a more 
sophisticated outer boundary treatment with a different damping profile can be found that will produce a 
numerical domain disturbance that can be controlled and made as small as desired. That is another study. 
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6.2 PML B on Four Sides 


In this section we consider simulations with PML B on all four sides, with all four corners. We consider only 
the case (M x , M y ) — (0, 0) with no convection. The complete domain is 

n c = [Xt - Du X r + D r ) x [Y b - D bl Y t + D t \, (69) 

with Xi = X r = Y b = Y t = 2, so that the inner numerical domain is the tiny 

n N = [-2,2] x [-2,2], (70) 

and the PML widths are all D; = D r = D b = D t = 9. The complete domain ilc is the same as the domain 
with BED outer boundaries that was used for both the single boundary computations reported in Table 1, 
and the blended corner computations reported in Table 2. The outer boundaries are at \x\ = 11 = \y\, so 
that the radiated signal will not reach the outer boundary by t = 10, and reflections from the outer boundary 
will return to the PML interfaces by t = 20. 


Table 4: Maximum absolute disturbance with PML on all four sides. 


PML B with Aa: = 1/24 = A y. At = 1/48. 


c 

t = 10 

t = 25 

t = 50 

0.001 

1.0653 x 10“ t2 

9.4618 x 10" 7 

3.4832 _y 

0.010 

1.0829 x 10” 11 

9.4618 x 10" 7 

3.4829“ y 

0.100 

1.2589 x 10~ 10 

9.4601 x 10- 7 

3.4806 _B 

1.000 

3.0775 x 10" y 

9.2995 x 10- 7 

3.5239 _B 

10.00 

2.0769 x 10 -Y 

2.8321 x 10" 7 

3.7187 -7 


The experiments reported in Table 4 handle the outer boundary by setting all of the data for all three 
primitive variables equal to 0 at \x\ = 11 = \y\, or with what we have called OBi in the previous section. We 
note immediately that at t = 10, before the signal reaches the outer boundary, the disturbance introduced 
by all four PMLs and corners is of the same order as the error at the same time with just one PML on 
the right, and even a bit smaller in most cases. At t = 25 the original signal has propagated through the 
PML to the outer boundary, and its reflection has returned all the way back across both the PML and the 
numerical domain. By t = 25 most of the accuracy of the PML has been lost due to this reflection from the 
outer boundary, as we have seen in the previous section. Continuing on to t = 50, the error grows another 
order of magnitude except for the heavily damped case with c = 10. With c = 10 the propagating signal is 
approximately O[10 -15 ] at x — xr = 4.75 in the PML B layer, and approximately O[10 -35 ] at x—xr = 9.75, 
with the signal being damped by approximately five order of magnitude over a PML length of about one and 
a quarter spatial units. With c = 1 the propagating signal is approximately O[10 -5 ] at x — xr = 4.75 in the 
PML B layer, and approximately O[10 - '] at x — xr = 9.5, with the signal being damped by approximately 
one order of magnitude over a PML length of about two and a half spatial units. For c = 10, with no flow, 
a PML B width of D = 9 appears sufficient to keep the numerical disturbance caused by the PML system 
from growing in time. For c = 1, with no flow, a PML width of from D = 27 to D = 30 would be required 
to reduce the propagating signal to O[10 -15 ]. Prohibitvely wide PML B layers would be needed to maintain 
any accuracy with smaller damping coefficients if such a crude outer boundary is used. We note also that the 
results reported in Table 4 show for PML B that with sufficient damping for the layer width, the disturbance 
appears stable in time. 


7 PML C with a Gaussian Jet Flow 

This series of computations is with a parallel Gaussian jet flow for the base solution, with 

Pb = 1, u b = 0.8 exp(— 36y 2 ), and; v b = 0. (71) 
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Note that the convection velocity is not uniform, though it is a parallel flow, and it does satisfy the Euler 
equations with uniform density. The maximum convection velocity is +0.8 at the center of the jet. In this 
case, the perturbation solution of the Linearized Euler Equations satisfies a more general equation form than 
with a uniform base solution, which produces the constant coefficient LEE. The same pressure source as 
before is used in the numerical domain but not in the PML domain. 


Table 5: Maximum absolute disturbance with a parallel jet. 
PML C with Ax = 1/16 = Ay, At = 1/32. 


c 

E C ( 4) 

E C { 8) 

0.001 

1.16155 x 10" y 

6.74813 x 10" y 

0.010 

1.16102 x 10" 8 

6.73600 x 10 -s 

0.100 

1.15577 x 10-' 

6.61588 x 10- y 

1.000 

1.10454 x 10" 6 

5.60609 x 10 _B 

10.000 

7.97674 x 10 _B 

3.06075 x 10 -b 

100.000 

3.92057 x 10 _t> 

6.58872 x lO” 5 


The simple PML C damping layer is used for the series of numerical experiments reported in Table 5. 
For these experiments we have included results with damping coefficient c = 100. Note that a coarser grid 
resolution is used for these experiments, with Ax = 1/16 = Ay and At = 1/32. The data in the first column 
is for t = 4, with numerical domain $l N = [—5, 8] x [—5, 5], and with PML C layer fi# = [8, 14] x [—5, 5], with 
interface X r = 8. The data in the second column is for t = 8, with numerical domain f I jv = [—9, 16] x [—9, 9], 
and with PML C layer il p> = [16, 30] x [—9, 9], with interface X r = 16. In both cases, the domains are large 
enough so that the propagated solution does not reach any outer boundary, either for the numerical domain 
or for the PML C layer. The data reported in Table 5 is for the maximum absolute disturbance introduced 
into the numerical solution by the PML C as measured at the interface, where the disturbance is greatest. 
As before, the data in Table 5 should be scaled by the pressure solution, and along the PML interface there 
are local wave crests with absolute magnitudes that are O[10 -4 ]. With a parallel jet steady mean flow, the 
disturbance of the numerical solution that is produced by introducing a PML C on the outflow boundary 
is comparable to the experiments with a uniform mean flow, ranging in absolute terms from O[10~ 9 ] to 
O[10 -5 ], and in relative terms from 0[1O -5 ] to O[10 -1 ]. As in the case of a uniform base flow, the disturbing 
effect of the PML is controlled by the size of the damping coefficient c, and the effect of the PML can be 
made as small as desired by taking c sufficiently small. For all but the largest two values of c, the effect of the 
PML is invisible at the scale of the solution. We see from these experiments that a PML can be used with 
nonuniform base flows, and the disturbance that is produced in the numerical solution can be controlled. 


8 Conclusions 

This paper presents the results of a series of numerical experiments to asses and compare PML boundary 
treatments. Three PML formulations have been used: PML A is an analytic treatment by Hagstrom [13]; 
PML B is from a characteristic analysis similar to that by Kami [16]; PML C is a familiar damping or 
sponge zone like that in [20] . PML A is formulated specifically for uniform base flows, and is somewhat more 
accurate in this case. PML B is simple, and overall is as effective as PML A, while it can be extended to 
variable coefficient and nonlinear systems. PML C is simple, robust, and the most easily generalized, but 
the disturbance that it produces in the numerical domain are a couple of orders of magnitude greater than 
from PML A or B. With a suitable layer width, any of these PML treatments can keep the errors that it 
introduces to less than 1% of the solution. All of the reported data is from computations with a uniform 
flow, except for one series of experiments with a parallel jet base flow. The damping profile a that we have 
used is smooth at the PML interface, and asymptotes to a = 1. The results from the reported numerical 
experiments show that for the Linearized Euler Equations: 

1. The level of distortion in the numerical solution that is introduced by each of the PML treatments is 
essentially linearly dependant upon the scaling coefficient c. The distortion levels introduced by the 
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PML treatments range from O[10 -10 ] to O[10 -1 ] in relative terms, and are controlled by the size of 
c. In this sense the PML treatments work well. PML A and PML B are similar in performance, and 
both produce disturbances that are orders of magnitude smaller than those produced by PML C. 

2. The blending of two PML treatments in a common corner can be done without introducing any 
additional distortion into the numerical solution domain. This can be done for a wide variety of PML 
treatments, including PML B and PML C, but it is not clear how this blending method would work 
for all PML treatments. 

3. The outer boundary of a PML layer can cause significant deterioration of the PML performance. Data 
was given only for PML B, which both damps to the right and amplifies to the left. Three outer 
boundary treatments were compared: OB\ which sets all data to 0 at the outer boundary; OB 2 which 
extrapolates R = P+u and v, and sets L = p — u = 0; and OB 3 which propagates R and v with interior 
differencing while setting L = 0. The errors introduced by the outer boundary conditions decreased 
linearly with the layer width, but they overwhelmed the PML treatments with the layer thicknesses 
that we considered except for the strongest damping with c = 10. 

4. The performance of PML C is essentially the same with uniform and nonuniform parallel base flows. 


The PML treatments that we have considered result either from operator approximation, which is difficult 
in the variable coefficient or nonlinear cases, or from eigenvector analysis, which is inherently one dimensional 
and misses the essentially multidimensional dynamic of systems like the Linearized Euler Equations in two 
or three space dimensions. All of the results that we have reported are dependant upon the accuracy of the 
code implementations that have been used, and while validation has been attempted, it is not a proof of 
coding correctness. Areas that need further investigation include NonReffecting Boundary Conditions for 
PML outer boundaries, or modifications to PML treatments that can control the error reflected from an 
outer boundary. 
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